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Abstract. Fully numerical mesh solutions of 2D and 3D quantum equations of 
£L|| Schrodinger and Hartree-Fock type allow us to work with wavefunctions which pos- 

sess a very flexible geometry. This flexibility is especially important for calculations 
' of atoms and molecules in strong external fields where neither the external field nor 

. the internal interactions can be considered as a perturbation. In the framework of this 

method we present various approaches to calculations of quasi-steady states of these 
systems in strong electric fields. These approaches are aimed at obtaining precise com- 
plex wavefunctions and corresponding eigenvalues in the form E = E$ — iT/2, where 
Eo is the real part of the energy and the value T/2 determines the lifetime of the state 
in relation to escape of electrons from the system. The applications for single-electron 
systems include correct solutions of the Schrodinger equation for the ion (energies 
and decay rates) and the hydrogen atom in strong parallel electric and magnetic fields. 
^ . Some results for the helium atom in strong electric fields are announced. 
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F^: I INTRODUCTION 

> ■ 

During the latter decade the interest to theoretical studies of atoms and molecules 
in strong external fields was strongly motivated by experiments with intense laser 
beams (electromagnetic fields with dominating electric component) and astronom- 
ical observations of white dwarfs and neutron stars (magnetic fields). The exper- 
imental availability of extremely strong electric fields in laser beams makes the 
theoretical study of various atomic and molecular species under such conditions 
very desirable. The properties of atomic and molecular systems in strong fields 
undergo dramatic changes in comparison with the field-free case. These changes 
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are associated with the strong distortions of the spatial distributions of the elec- 
tronic density and correspondingly the geometry of the electronic wavef unctions. 
This complex geometry is difficult for its description by means of traditional sets 
of basis functions and requires more flexible approaches which can, in particular, 
be provided by multi-dimensional mesh finite-difference methods. 

Many results of the experiments with intense laser beams can be considered from 
the point of view of the behaviour of atoms and molecules in intense static electric 
fields, especially when the frequency of the radiation is low. (The low-frequency 
behaviour and the limits of this region are analysed by [1]). 

The most advanced theoretical studies of effects in strong electric fields were tra- 
ditionally concentrated on a hydrogen atom. Other atomic and molecular systems 
in strong electric fields are much less studied. Many sophisticated theoretical meth- 
ods developed for the hydrogen atom cannot be simply applied for other atoms and 
molecules. This circumstance is an argument for development of more universal 
theoretical and computational methods for atoms and molecules in strong electric 
fields. 

Quasi-steady states of hydrogen atoms in strong electric fields were studied pre- 
cisely in many theoretical works (see [2-11] and references therein). Some of 
these works are based on separation of variables in parabolic coordinates in the 
Schrodinger equation for the hydrogen atom. This separation of spatial variables 
for atomic and molecular systems in external electric fields is possible only for one- 
electron atoms with the Coulomb electron-nucleus interaction. Non-hydrogenic 
systems or systems with non-Coulomb interaction do not allow this separation 
of variables. Their theoretical studies require a developed technique of solving the 
Schrodinger and similar equations with non-separable variables. One of the possible 
ways for solution of this problem consists in the application of mesh computational 
methods for solving these equations. The mathematical problem consists here in so- 
lution of partial differential equations for systems with discrete quasi-steady states 
lying on the background of the continuous spectrum. 

In the first part of this work we present applications of our mesh method for 
solving Schrodinger equations with non-separable variables for quasi-steady states. 
The method is based on the technique of a mesh numerical solution of Schrodinger 
and Hartree-Fock equations with non-separable variables for steady states [12-17]. 
The most of the applications for the discrete states was focused on atoms in strong 
magnetic fields [14-20]. In this paper we present computational approaches which 
we have developed for quasi-steady states in external electric fields [21-23] and 
apply them to several single electron systems. At the end of the paper we announce 
some of results of our current work on the helium atom in strong electric fields. 



II FORMULATION OF THE PROBLEM AND THE 2D 
MESH COMPUTATIONAL METHOD FOR 
STATIONARY STATES 



We carry out our mesh solution in the cylindrical coordinate system (p, <p, z) with 
the z-axis oriented along the electric field. After separation of the coordinate the 
Hamiltonians of the single-electron problems considered below take the form 
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- the molecular ion in electric and magnetic fields parallel to the molecular axis. 
Here and in the following the atomic units m e — h — e — 1 will be used, including 
the magnetic field strength 7 measured in units B = hc/eal = 2.3505- 10 5 T, 7 = 
B/Bq. F is the electric field strength multiplied by the charge of the electron. 
F — 1 corresponds to 51.422V/A. The value m is the magnetic quantum number 
and s z = ±i is the spin ^-projection. R is the internuclear distance for Hj" molecule. 
The Hamiltonian (2) does not include the internuclear repulsion and its eigenvalues 
are the electronic energies E e which are connected with the total energy of H<j~ by 
the formula E = E e + 1/R. 

Our two-dimensional mesh computational method which we apply for obtain- 
ing eigenvalues of the Hamiltonians (1) and (2) in the form sufficient for case of 
stationary states (i.e F — 0) of single-electron systems is described in refs. [13,14]. 
For these systems highly precisely solutions can be obtained by solving Schrodinger 
equations in finite spatial domains Vt (with simple boundary conditions V'lan = 
or dip/dn\ dn = 0) with negligible errors for domains of reasonable dimensions (see 
analytical estimations of errors in ref. [14] and results of numerical experiments 
in ref. [13]). When employing the Richardson's extrapolation of the energy values 
to the infinitely small mesh step h — > very precise results can be obtained on 
uniform meshes with relatively small numbers of nodes [13,14]. Two important 
problems have been solved for the following development of the multi-dimensional 
Hartree-Fock method for many-electron systems: 1. An enhancement of the pre- 
cision of single-electron wavefunctions due to their more complicated geometry in 
comparison with the one-electron case, and 2. Obtaining correct mesh represen- 
tations of Coulomb and exchange potentials. The first of these problems can be 
solved more or less simply by means of non-uniform meshes with the distribution of 



nodes concentrated near the nuclei. The second and most complicated problem was 
initially solved in our first works on the helium atom in magnetic fields [15,16] by 
means of a direct summation over the mesh nodes. But this direct method is very 
expensive with respect to the computing time and due to this reason we obtained 
in the following works [17-20] these potentials as solutions of the corresponding 
Poisson equations. The problem of the boundary conditions for a Poisson equation 
as well as the problem of simultaneously solving Poisson equations on the same 
meshes with Schrodinger-like equations for the wave functions ip^z.p) has been 
discussed in ref. [16]. 

The simultaneous solution of the Poisson equations for the Coulomb and ex- 
change potentials and the Schrodinger-like equations for the wave functions ip^z, p) 
is a complicated problem due to a different asymptotic behaviour of the wavefunc- 
tions and potentials. The wavefunctions of the bound electrons decrease exponen- 
tially as r — > oo (r is the distance from the origin). This simplifies the problem 
of the solution of equations for wavefunctions in the infinite space because it is 
possible either to solve these equations in a finite domain Q (as described above) 
or otherwise to solve these equations in the infinite space on meshes with expo- 
nentially growing distances between nodes as r — > oo. On the contrary, solutions 
of Poisson equations for non-zero sums of charges decrease as 1/r as r — > oo. In 
result, every spatial restriction of the domain fl introduces a significant error into 
the final solution. In our approach we address the above problems by using special 
forms of non- uniform meshes [17]. Solutions to the Poisson equation on separate 
meshes contain some errors 5p associated with an inaccurate description of the 
potential far from the nucleus. _ However, due to the special form of the function 
Sp(h) for these meshes (where h is a formal mesh step) the errors do not show up 
in the final results for the energy and other physical quantities, which we obtain 
by means of the Richardson extrapolation procedure (polynomial extrapolation to 
h = [13,14]). The main requirement for these meshes is a polynomial increase 
of the actual mesh step h(r) when r — > oo. Moreover, this behaviour can be only 
linear one, i.e. h~ l (r) = 0(l/r) as r — > oo. The error of the mesh solution in this 
case has the form of a polynomial of the formal step of the mesh h — 1/N, where 
N is the number of nodes along one of the coordinates. In practical calculations 
these meshes are introduced by means of an orthogonal coordinate transformation 
from the physical coordinates x p to the mathematical ones x m made separately for 
p and z. The numerical solution is, in fact, carried out on uniform meshes in the 
mathematical coordinates x m . The characteristic feature of these meshes consists 
of rapidly increasing coordinates of several outermost nodes when increasing the 
total number of nodes and decreasing the actual mesh step in the vicinity of the 
origin. 

The methods described above for the two-dimensional case can be applied also 
to three-dimensional Schrodinger and Hartree-Fock problems [24,25]. 



Ill QUASISTATIONARY STATES 



The most common approach to the mathematical description of the quasista- 
tionary states of quantum systems consists in employment of the complex form of 
the energy eigenvalues 

E = E - iT/2 

where the real part of the energy Eq is the centre of the band corresponding to 
a quasistationary state and the imaginary part T/2 is the half- width of the band 
which determines the lifetime of the state. In this approach one may consider 
quasistationary states of quantum systems similarly to the stationary ones. For 
systems (1) and (2) at F > the electron can leave the system in the direction of 
positive z and the behaviour of the wavefunction on this semi-axis determines the 
main features of the behaviour in the external electric field F. The mathematical 
problem consists in the solution of the Schrodinger equation 

Hif) = Eif) (3) 

at F > for resonance states which are inheritors of some discrete (at F = 0) states. 
The wavefunctions of these states must describe the process of separation of the 
electron from the system. These wavefunctions can be distinguished either by the 
explicit establishing the boundary condition of the outgoing wave or by means of a 
complex coordinate transformation. The latter transforms the oscillating outgoing 
wave into an exponentially decreasing wavefunction for which the simple Dirichlet 
boundary condition if) — > as z — ► +oo can be established. 

From the mathematical point of view the problem under consideration consists in 
obtaining solutions of the single-particle Schrodinger equation for this electron with 
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FIGURE 1. The potential energy for the hydrogen atom V(p — 0, z) in the external uniform 
electric field. 



the correct asymptotic behaviour of the wavefunction as an outgoing wave. Cur- 
rently we have three different possibilities for obtaining solutions with the outgoing 
wave asymptotics realised in our computational program: 

1. Complex boundary condition method. This method is described in detail in 
ref. [16,22]. The method is based on the idea that the single-electron Schrodinger 
equation for a finite system can be solved with the arbitrary precision in a finite 
area both for stationary and for quasi- stationary eigenstates. The case of stationary 
states is considered in [13,14]. We discuss the approach for the quasi- stationary 
states following ref. [22]. Figure 1 presents the potential curve for the simplest 
Hamiltonian (1). Analogously to [14,27] the calculations can be carried out in 
an area Q which is finite along the direction z. For the coordinate z we have 
used uniform meshes. The boundary of the -L 2 _ for z < (F > 0) 

(Figure 1) is determined from the condition of small values of the wavefunction on 
the boundary and, therefore, small perturbations introduced by the corresponding 
boundary condition [14]. The values of the wavefunction on the opposite boundary 
of the area (z = L z+ ) cannot be excluded from the consideration. We consider 
non- stationary states corresponding to the process when an electron leaves the 
system in the direction z — > +oo. Thus, an outgoing wave boundary condition 
is to be established on z = L z+ . The form of this boundary condition can be 
derived from the asymptotic behaviour of the wavefunction for z — > +oo. In this 
limit the asymptotics of the real system can be replaced by the asymptotics of 
the wavefunction of a free electron in the uniform electric field. Solutions of the 
Schrodinger equation for this system can be written as 

Mz) = AM(0e~^\ + ( 2F ) 1/3 ( 4 ) 

where M(£) and @(£) are the modulus and phase of the Airy function, A is a 
constant [26]. The asymptotics of ip for z — > +oo can be obtained from eqs. 
(10.4.78) and (10.4.79) of [26] in the form 

m = -^r 1/4 ex P + r\e /2 ) + o (r 13/4 ) (5) 

Taking into account that d£/dz = (2F) 1 / 3 and ^ 1 / 2 (2F) 1 / 3 = [2{E + Fz)} 1 / 2 we 
have from (5) 

where k — [2(E + Fz)} 1 / 2 is the wavenumber. Equation (6) allows us to establish 
the following outgoing wave boundary condition on the upper (in the z direction) 
edge of the region 



This approximate boundary condition is derived from the asymptotics of the wave- 
function of a free electron in a uniform electric field and in the limit L z+ — ► +00 can 
be considered as exact one. On the other hand, our numerical experiments show, 
that errors caused by the approximate nature of this asymptotics are important 
only for very short regions Q. Solving the Schrodinger equation with the bound- 
ary condition (7) established on a reasonable distances L z+ from the origin of the 
system one obtains precise complex eigenvalues of the energy and corresponding 
wavef unctions. 

This straightforward approach enables obtaining precise results both for atoms 
and molecules from weak to moderate strong fields (for instance for the ground 
state of the hydrogen atom up to F = 0.20 — 0.25 a.u.). For stronger fields the 
precision of this method is limited by the precision of the mesh representation of 
the boundary condition (7). 

2. Classical complex rotation of the coordinate z in the form z — > ze tB . Oppo- 
site to the boundary condition approach, this method does not require establishing 
boundary conditions dependent on the energy and can be used with the zero bound- 
ary conditions for the wavefunction. This simplification is especially important for 
applications to multi-electron systems. In this approach we have obtained precise 
results for atomic systems in strong fields from the lower bound of the over-barrier 
regime up to super-strong fields corresponding to regime |Re.E| << |Imi£| [23]. For 
atoms in weak fields the applicability of the method is limited by numerical errors 
in the imaginary part of the energy. On the other hand, this method cannot be 
immediately applied to molecular systems in our direct mesh approach [23,28]. 

3. Exterior complex transformation of the coordinate z. The exterior complex 
scaling [28] combines many advantages of the boundary condition and complex 
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FIGURE 2. The integration paths for the complex coordinate z. Dashed line is the path at the 
traditional rotation z — > ze tn . Solid line is the path corresponded to the exterior curved path 
coordinate transformation. 



rotation method. In its initial form [28] it consists in the complex rotation of a 
coordinate e.g. r around a point r lying in the exterior part of the system, i.e. 

t — ► i"o + ( r — ^o)e* for r > r (8) 

and leaves intact the Hamiltonian in the internal part of the system. The latter 
circumstance allows us to employ this transformation both for atoms and molecules 
in arbitrary electric fields. In addition this transformation does not contain energy- 
dependent boundary conditions as well as the classical complex rotation. On 
the other hand, the exterior complex scaling being introduced into our numeri- 
cal scheme in the form (8) leads to the same numerical problems at very strong 
fields as well as the boundary condition method due to the nonanalytic behaviour 
of this transformation at the point r (or z in our case). 

In our numerical approach we solve the latter problem by means of a transfor- 
mation of the real coordinate z into a smooth curved path in the complex plane z 
(Figure 2, see details in [23]). This transformation leaves intact the Hamiltonian 
in the internal part of the system, but supplies the complex rotation of z (and the 
possibility to use the zero asymptotic boundary conditions for the wavefunction) in 
the external part of the system without any loss of analycity. The transformation 
can be applied both for atoms and molecules and provides precise results for fields 
from weak up to super-strong including the regime |Re.E| << |IhlE| [23]. 

The three methods presented above have different but overlapping regions of 
their most effective application. This allows their combined using for the recipro- 
cal control in applications. The numerical results obtained by all three methods 
coincide in the limits of their applicability and are in agreement with numerous 
published data on the hydrogen atom in electric fields (see e.g. [5-8]). 

IV SELECTED PHYSICAL RESULTS 

Using the boundary condition approach enables one to obtain both the values 
of the energy and other observables and the wavefunctions in their mesh represen- 
tation. An example of such a wavefunction of the ground state of the hydrogen 
molecular ion is presented in Figure 3. The electron energy E c obtained as a so- 
lution of (3), (2) allows one to determine potential curves E(R) for the molecule 
as a whole by using the formula E(R) = E c + 1/R. These potential curves and 
corresponding values of T/2 are presented in Figure 4. One can see in Figure 4 
that when growing electric field strength the minimum on the curve E(R) shifts 
to the right, becomes more shallow and flat, and at F c pa 0.065 the minimum dis- 
appears. The dependence of the location Ro of the minimum on these curves with 
corresponding values of the energy E and half-widths of the level are presented in 
Table 1. (For comparison, at F = published results for the equilibrium point are 
Rq = 1.997193 and E = —0.6026346 [29]. For the potential curve as a whole our 
numerical values at F = coincide with data by [30].) One can see in Table 1 



and Figure 4 that for the ground state of the hydrogen molecular ion at equilib- 
rium internuclear distances there is a marked probability of its decay through the 
separation of the electron. Analogous calculations for two excited states 2pir n (the 
lowest state with \m\ = 1) and 3d5 g (the lowest state with \m\ = 2) [22] show that 
for the field strengths at which their potential curves have minima this process has 
a low probability. 

The critical value of the maximal electric field F c = 0.065a.u. = 3.3V/A when 
the molecule Hg" can exist is in a good agreement with experimental results by 
Bucksbaum et al [31]. According this work molecule may exist in laser beam 
fields with intensity less than 10 14 W/cm 2 which corresponds to 3V/A, and does 
not exist in more intensive fields. 

Other theoretical explanations of the rupture of molecule in intensive laser 




FIGURE 3. Real and imaginary parts of the wavefunction of the H2 molecule in a longitudinal 
electric field F (a.u.). 




FIGURE 4. Left - Potential curves for the ground state of the Hj" molecule in the longitudinal 
electric field F (a.u.). Right - Half-widths of the ground state energy level of the molecule Hj" in 
the longitudinal electric field F as a function of the internuclear distance R (a.u.). 



TABLE 1. Equilibrium internuclear distances and corresponding en- 
ergies and half- widths of the energy of the ground state of H J molecule 



in longitudinal electric fields. 
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0.00 


1.997 


-0.60264 


0.04 


2.062 


-0.60686 


1.23(-19) 


0.01 


2.001 


-0.60289 


0.05 


2.112 


-0.60943 


7.32(-15) 


0.02 


2.012 


-0.60366 


0.06 


2.198 


-0.61285 


1.49(-11) 


0.03 


2.031 


-0.60497 


0.065 


2.28 


-0.61501 


3.94(-10) 



beam fields in works [31-33] were based on conception of deformation of potential 
curves which results in coupling of the ground state lscr g with 2pa u state [32]. This 
results in a possibility of rupture of the molecule after absorption some photons. 
On the other hand, the results obtained above enable us to analyse this process 
from simpler point of view. When the frequency of the radiation is low enough the 
consideration of the process as extension and rupture of the molecule in a strong 
static electric field is a quite adequate model [1,32]. The conception of coupling 
states might be substituted by an exact numerical calculation of one (as it was 
done above) or several (if several states are considered) dependencies E(R). This 
numerical calculation is equivalent to a traditional calculation of the energy of the 
state when taking into account coupling with all the states of the same symmetry 
(corresponding to the symmetry of the initial physical problem). 

The calculations presented above were carried out for the ground state of 2D 
Hamiltonian. In fact, our method is not restricted by these states. Computa- 
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FIGURE 5. Left - The electron energies E e of 2sa g and 3da g states of the HJ molecular ion 
as functions of the internuclear distance R in the field-free space (broken), in the magnetic field 
7 = 0.02 parallel to the molecular axis (full) and in the electric field F = 0.002 parallel to 
the molecular axis (dotted). Right - E e (F) for the doublet of 2sa g and 3da g states of at 
R = 4.0536a.u. for 7 = (broken) and 7 = 0.02 (full). 



TABLE 2. The ground state of the hydrogen atom in parallel electric and magnetic fields 
at F = 0.1, 1, 5. 
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-7.82561 


1.04701(-2) 



tional algorithms of our program of atomic and molecular mesh calculations AT- 
MOLMESH are constructed so that calculations are being carried out for a state 
with prescribed spatial symmetry and with the electron energy nearest to the given 
initial approximation E\>. Thus, there is no difference between calculations for the 
ground and for excited states and calculations can be fulfilled for various states of 
the same spatial symmetry including even degenerate states (when applying some 
special technical methods). For instance, potential curves 2sa g and 3dcr g of the 
field- free fi~2~ have exact crossing near R = 4.05 [30]. In magnetic fields parallel 
to the molecular axis this pair of states forms an avoided crossing. Our method 
permits calculations for both exact and avoided crossings at F = as well as at 
F^Oas presented in Figure 5. 

The second singe-electron system which we present here is the hydrogen atom in 
parallel electric and magnetic fields [23]. Some numbers obtained for this system 
are shown in Table 2. 
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FIGURE 6. Left - The total energies of 2 1 So state of the helium atom dependent on the electric 
field strength. Right - Half-widths of this energy level. 



Concluding this section we announce some preliminary results on the helium 
atom in strong electric fields. These results are obtained in a multi-configurational 
(CI) mesh calculation. The calculations for separate meshes are analogous to de- 
scribed above and are carried out by means of the curved path exterior complex 
transformation. Our preliminary results for the ground state of the helium atom 
agree with data from recent works [34,35]. Our results for the excited state 2 1 S 
are shown in Figure 6. 

V CONCLUSIONS 

In this communication we have presented a 2D fully numerical mesh solution 
approach to calculations for atoms and simple diatomic molecules in strong ex- 
ternal electric and magnetic fields. For single-electron systems in external electric 
fields we can apply three different methods of calculation of the complex energy 
eigenvalues for atom-like systems and two methods for molecules. These methods 
have different but overlapping regions of their most effective application. This al- 
lows their combined using for the reciprocal control in applications. The complex 
boundary condition method is the most straightforward and is the most reliable 
in this sense. For relatively weak fields it allows obtaining results on the meshes 
with the lesser number of nodes than the curved path complex transformation. On 
the other hand, the boundary condition method loses its precision and stability at 
extremely strong electric fields and, second, this method contain the energy of the 
electron in its formulation. The latter feature can be some obstacle in its applica- 
tion to problems more complicated than the single-electron Schrodinger equation. 
The curved path complex transformation method is the most general of the three 
and the most prospective for the following applications. The only shortage of this 
method is that this is less precise and less stable for the atom-like systems at 
extremely strong electric fields (|Imi?| >> |Re.E|) than the traditional complex 
rotation. Thus, the latter method is the most convenient for atom-like systems at 
such extremely strong fields. 

The mathematical technique developed for solving Schrodinger equations for 
quasi-steady states allowed us to obtain a series of results for the hydrogen atom in 
parallel electric and magnetic fields and for the ion in strong electric fields. The 
following applications of these methods are associated with the CI approach which 
is now in the process of development and testing on the problem of the helium 
atom in strong electric fields. We present preliminary results for the 2 1 S state of 
this atom. 
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